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The effects of mid-range repulsion in Lattice Boltzmann models on the coalescence/breakup be- 
haviour of single-component, non-ideal fluids are investigated. It is found that mid-range repulsive 
interactions allow the formation of spray-like, multi-droplet configurations, with droplet size di- 
rectly related to the strength of the repulsive interaction. The simulations show that just a tiny 
ten-percent of mid-range repulsive pseudo-energy can boost the surface/ volume ratio of the phase- 
separated fluid by nearly two orders of magnitude. Drawing upon a formal analogy with magnetic 
Ising systems, a pseudo-potential energy is defined, which is found to behave like a quasi-conserved 
quantity for most of the time-evolution. This offers a useful quantitative indicator of the stability 
of the various configurations, thus helping the task of their interpretation and classification. The 
present approach appears to be a promising tool for the computational modelling of complex flow 
phenomena, such as atomization, spray formation and micro-emulsions, break-up phenomena and 
possibly glassy-like systems as well. 

PACS numbers: 

I. INTRODUCTION 

In the last two decades, the Lattice-Boltzmann (LB) approach has emerged as powerful mesoscopic alternative to 
classical macroscopic methods for computational hydrodynamics [lM3|. The pseudopotential method, put forward 
a decade ago by X. Shan and H. Chen to endow Lattice Boltzmann models with potential energy interactions, 
is one of the most successful outgrowths of basic LB theory 0, [Bf. The Shan-Chen (SC) model is based on the 
idea of representing intermolecular interactions at the mesoscopic scale via a density-dependent nearest-neighbour 
pseudopotential ip(p). Despite its simplified character, the SC model provides the essential ingredients of non-ideal 
fluid behaviour, namely a non-ideal equation of state and surface tension effects at phase interfaces. Due to its 
remarkable computational simplicity, the SC method is being used for a wide and growing body of complex flows 
applications, such as multiphase flows in chemical, manufacturing and geophysical problems. 

To date, the overwhelming majority of Shan-Chen applications have been performed within the original formulation, 
whereby only first-neighbor attractive interactions are included. This entails a number of limitations, primarily the 
impossibility to tune the surface tension independently of the equation of state. This limitation has been recently 
lifted by introducing second- neighbor repulsive interactions [6| . Besides offering an independent handle on the surface 
tension, it has been observed that second-neighbor (mid-range) repulsion may disclose an entirely new set of physical 
regimes, particularly the onset of metastable multi-droplet configurations, which would be impossible to obtain with 
short-range attraction alone. These configurations result from the existence of energy barriers (mid-range repulsion) 
which slow- down/ arrest the dynamics of coarsening/phase-separation [7|, [8[ In this work, we provide a quantitative 
exploration of the basic mechanisms behind this physically enriched scenario. To this aim, we investigate the structural 
properties of multi-droplet configurations, as well as their energetics, as a function of the main parameters of the model, 
mainly the strength of the repulsive interactions. Upon progressive switching of this paramater, the system is found 
to move from a single-droplet phase-separated fluid, to a multi-droplet metastable configurations, all the way up to a 
quasi-ordered crystal-like structure. 
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II. STANDARD SHAN-CHEN MODEL 



The standard lattice Boltzmann (LB) equation with pseudopotential interaction can be expressed as follows, 

fi(x + Ci,t + 1) - fi{x,t) = -L0(fi - ff) + Fi(x) (1) 

where fa is the probability density function of finding a particle at site r at time t, moving along the i — th lattice 
direction defined by the discrete speeds q with i = 0, The left hand-side of (pQ) stands for molecular free- 

streaming, whereas the right-hand side represents the time relaxation (due to collisions) towards local Maxwellian 
equilibrium. Finally, Fi represents the total volumetric body force. In particular, we shall use a dynamic mean- field 
term connected with bulk particle-particle interactions. The macroscopic density p and velocity u are given by [3| 

b 

p{x 1 t) = Y J h (2) 



p{x , t)u{x , t) = ^ %fi ( 3 ) 

i=0 

The equilibrium distribution function is calculated in order to make the collision operator conserve mass and 
momentum: a common choice that satisfies the above constraints is the following 

^ = VHp(l + ^ • u + Jjft ■ uf - (4) 
The Fi term in ([1]) represents the phase interaction, 

b 

F = -Gi/j(x) ^2 w i^(% + (5) 

i=0 

in which ip(p) is the local pseudopotential governing the interaction and Wi are statistical weights which will be defined 
in the following. The expression of ijj{p) by Shan and Chen is the following, 



iP = ^(l-e-p/P°) (6) 

In this model, phase separation is achieved by imposing a short-range attraction between the light and dense phases. 
Indeed, such short-range attraction is responsible for the growth of density contrasts through a dynamical instability 
of the interface. In real fluids, such instability is tamed by hard-core repulsion, while in the SC model such hard-core 
repulsion is not included, for it would impose significant penalty on the time-marching procedure, and is replaced 
instead by a saturation of the attractive interactions above a given density threshold, po. Expanding i^, ([5]), in terms 
of (5, we find, to fourth order [9| 



F = -dc 2 s G^Vi) - C 2 c A s G^VV 2 i) (7) 

where C\ and C2 are lattice-specific numerical factors. The first term is responsible for the non-ideal part of the 
corresponding equation of state: 



P = pel + \cxc\Gp V> 2 . (8) 
The second term in eq. ([7j) is the inherent surface tension in the SC model which yields 

7 = -^-<4 J Jd v n 2 dy (9) 

Considering the interface-equilibrium condition, ^ = ^f^Jr = 0, we find the critical condition for phase separation, 
G < G cr — —4.0, p cr = po ln2 
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FIG. 1: (Color Online) Two-belt lattice for force evaluation. Each node is labelled by the corresponding energy |cfci| 2 . Belt 1 
contains eight speeds and two energy levels (1, 2). Belt 2 contains sixteen speeds, distributed over three energy levels (4, 5, 8) 

III. SHORT AND MID-RANGE INTERACTIONS 

Our model is based on the interaction between each particle and a set of 24 surrounding neighbours, distributed 
over two Brillouin zones (belts for simplicity). The interaction force in (pp) reads as follows 

= ^2 G k^k(p(x))^2pkAi^k(p(x^c ki )) (10) 

k=s,m i—1 

where the index k = s, m labels the short and mid range belts respectively, whereas c%i denotes the i-th set of 
discrete speeds belonging to the k-th. belt. The pseudo-potential force consists of two separate components F(x, t) = 
F s (x,t) + F m (x,t), defined as follows: 

bi 

F a (x,t) = G 1 *l>(x,t)^w i 'il>(x i ,t)ciAt (11) 

2=1 

bi 

+ G 2 ^{X, t) y^Paj^Xai, t)c si At 
i=l 
b 2 

F m (x,t) = G 2 1p(x, t) ^2 Pmi^(Xmi,t)c mi At 

2=1 

In the above, Wi are the weights of the first belt of neighbours, the same as in the standard SC model; the indices 
k = s,m refer to the first and second Brillouin belts in the lattice, and c^, Pki are the corresponding discrete speeds 
and associated weights, reported in Tab. 1. Finally Xk% = x + At are the displacements along the i-th direction in 
the k-th belt. 



E(8) 


Psi -- 


= P (1) = 


= 4/63 , 


i = l,4 


Psi -- 


= p(2) = 


= 4/135 , 


z = 5,8 


Pmi 


= p(4) 


= 1/180 , 


i = l,4 


Pmi 


= P(5) 


= 2/945 , 


i = 5,12 


Pmi 


= p(8) 


= 1/15120 , 


z = 13,16 



TABLE I: Links and weights of the two-belt, 24-speed lattice [To|, [Tl[. 

Note that G is a measure of potential to thermal energy ratio, and positive (negative) G correspond to repul- 
sion(attraction) respectively. The first belt is discretized with 9 speeds (b\ = 8), while the second with 16 (62 = 16) 



4 



and the weights are chosen in such a way as to fulfill the following normalizations [l(J [H| : 



b 2 



^2 Wi = 5Z Psi + ^Vmi = 1 (12) 



z=0 



b\ b\ 62 

^Psicti + ^2PmiC 2 mi = C 2 S (13) 



i=l i=l i=l 



where c 2 s = 1/3 is the lattice sound speed. Note that the present set of discrete speeds and weights secures 8 — th 
order isotropy in the force evaluation. The pseudo-potential ip(x) is taken in the form first suggested by Shan and 
Chen [4|], ip[p] = ^/po(l — e~ p / p0 ) where po marks the density value (critical) at which non ideal-effects come into play 
and it is fixed to p Q = 1 in lattice units. Taylor expansion of (fTTj) to second-order delivers the following non-ideal 
equation of state (EOS) 



P^P/cl = P + { -^±^-^(x,t) (14) 

where gu = Gk/c 2 s are normalized coupling strengths. Further expansion of eq. ([TT]) to fourth-order provides the 
following expression for the surface tension 



(Gi + f G 2 



/CO 
\d y ^)\ 2 dy (15) 
-co 



where y runs across the phase interface. This is the analogue of eq. (j9j), with the correspondence : C\G (Gi + G2) 
andC 2 G^(Gi + ^G 2 ). 



IV. NUMERICAL RESULTS 



With two parameters at our disposal, G\ and G2, the present model allows a separate control of the equation of 
state and surface tension, respectively. In particular, as shown in previous work [12|, the non-ideal part of the equation 
of state depends only on A\ = G\ + G2, whereas surface tension effects are controlled by the combination G\ + ^G^. 
Since in the vicinity of 7 — » higher order terms come into play, it proves expedient to define a new coefficient 

A 2 = Gi + AG 2 (16) 

where the numerical factor A plays the role of a renormalisation parameter, whose departure from zeroth-order value 
^ is a measure of the influence of the higher-order terms. Comparison with numerical results shows that A w 3/2 
provides satisfactory agreement, see eq. (fT8]) . This shows that, at a given value of A\ (i.e. given density ratio between 
the light and dense phase), mid-range repulsion (G2 > 0) is expected to lower the surface tension of the fluid, thereby 
facilitating the formation of multi-droplet configurations with higher surface/ volume ratioes than the standard Shan- 
Chen model. Thus, the mid-range potential is expected to act as a "surfactant" fl3l~Tl5| . where "surfactant" indicate 
that true surfactants can be transported by many different mechanisms, locally changing the surface tension of the 
fluid, which is not what the mid-range repulsion in the present work does. Numerically, the role of the mid-range is 
to add higher-order derivative to standard interaction force, which provides more isotropy and enables control of the 
equilibrium surface tension. In order to explore this scenario, we have simulated droplet formation by integrating the 
LBE Eq.(pQ) in a 2D lattice using the nine- speed 2DQ9 model US-El, out of a noisy density background (Sp/p ~ 0.01) 
with initial density pi n = poln2 + Sp in a periodic domain. In all simulations, r = 1. We have performed a systematic 
scan over the force strength, by changing G\ and G2 so as to keep A\ = —4.9 while increasing A2 above the Shan- 
Chen value A2 = A\ = —4.9. All simulations have been performed with a resolution of 512 2 grid points, and a total 
simulation time t = 500000. 

In fig. [2j some snapshots of the density at final time are shown for different value of A2. In fig. [3j it is shown 
the number of droplets, at the end of the simulation, as a function of A2. The simulations show a threshold in 
phase-separation as G2 increases towards a critical value, A2 C ' beyond ^2 C , the density field exhibits numerous stable 
droplets, distributed according to a quasi-ordered configuration, somehow reminiscent of a crystal-like configuration 
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FIG. 2: (Color online) Spatial distribution of the fluid density. The formation of a large number of droplets with increasing A 2 
is well visible, [a]: Standard Shan-Chen, A x = -4.9, A 2 = -4.9 nx = ny = 512, t = 500000; [b]: Multi-droplet, A 2 = -2.85, 
nx = ny = 512, t = 500000; [c]: Multi-droplet, A 2 = -0.8, nx = ny = 512, t = 500000; [d]: Multi-droplet, A 2 = 0.15, 
nx — ny — 512, t = 500000. The right panel shows the Fourier spectrum of density fluctuations. Such spectrum, initially a 
white noise, evolves towards a shape peaked at the (inverse) size of the droplets. These Fourier spectra show that small-scale 
contribution is significantly higher when increasing the mid-range repulsion, that is A 2 , indicating the formation of long-lived 
metastable states in the form of small droplets. In particular in the last picture (d), at the end of the simulation there is a 
clear peak at R ~ L/2k w 30. 



with defects. The numerical value of A^ c can be roughly estimated by noting that, to fourth-order in the lattice 
spacing, the total force due to intermolecular interactions, F tot — F s + F m , is given by: 

Flat = - (ciA^Vip + ^ VVA^ (17) 

where, as previously mentioned, A\ = G\ + G2 controls the magnitude of the phase separation (liquid to gas density) 

2 

and A2 = G\ + AG2 is directly linked to the surface tension. A dimensional argument gives I 2 — ^fj^ — thus 
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FIG. 3: Number of droplets versus A2 after a short time t = 10000 and final time t = 5 X 10 5 . The vertical line denotes the 
transition zone from the Multi-Droplet to the Emulsion region, in correspondence with the theoretical spinodal point A2C = 0. 
The standard SC single-droplet region is associated with A2 — >> —4.9 (G2 —> 0). For < G2 < G2C metastable multi-droplet 
configurations are found, which tend, nevertheless, to the single- droplet equilibrium configuration after a sufficient long time. 
For A2 > 0, the relaxation time associated to the decay to this equilibrium state becomes formally infinite (no changes in time 
for all observables) , indicating that the non-equilibrium phenomena that sustain these metastable states experience very slow 
dynamics. The solid lines represent two different fits for the two regions, the multi-droplet and the emulsion one. Respectively, 
they are given by: n(t) = (1 - £)~ p(t \ with P(t) = i+t/iot eap and f = G 2 /G 2c ] n(f) = af with a = 10000. 



yielding 

1 AXT / / 

(18) 




where A = 3/2 has been used in the rightmost expression. In the above, I is the typical size of a nucleus and l\ is the 
typical single-droplet size for the Shan-Chen case. With this choice of A, the resulting spinodal value, at which I — >> 0, 
turns out to be G2 C = 9.8, corresponding to A2 C = 0.0. For this value, the coefficient in front of the second term in 
Eq. (fT7|) vanishes, thus signaling the onset of a phase-transition. This value is found to be in good agreement with 
the numerical simulations, which indicate complete nucleation starting around a value of A2 ~ 0, as shown in Fig. 
[U The number of droplets in the first region, called Multi-droplet region, can be described as a function of time as 
n (t) ~ (£ — l)~ p ^\ with p(t) = 1+t /i Qt — , where £ = G2/G2C and t cap — This relation can be related to simple 
statistical physics arguments [12j. The region after the transition, where nucleation takes place, has been fitted by 
a simple linear function n(£) = a£, where a = 10 4 . It is worth mentioning that the same linear behaviour in the 
emulsion region is also observed for a coarser domain [12|. However, the coefficient a is not universal, as it depends 
on the domain size. This may be related to the breakdown of scale invariance of phase-separating fluids as observed 
in Hp. 

It is instructive to inspect the spatial distribution of the phase-separated fluid as the A2 parameter is increased. 
In the present model, as well as in the standard Shan-Chen, phase separation starts immediately and spontaneously, 
once the parameters are chosen in the critical range: in the Shan Chen model and in the two-belt (with A2 below the 
critical value), small droplets coalesce in larger droplets of increasing size, until only a few of them, or even just one, 
are left. This is the spatial configuration which minimizes the surface energy expenditure. To study this spontaneous 
coalescence and its relation to the model parameters, a Fourier analysis of the density field has been conducted, based 
on the structure factor S(k,t): 



[p{x,t)-p{t)]e' 



k x 



(19) 



where k = (^)(/, m), x is the lattice point, L is the linear lattice size (= 512 in our case), N = L 2 is the total number 
of grid points; p(x, t) is the density field at time t and ~p(t) is the average density field at time t. It is possible to 

average the structure factor in k space, as follows: S(k,t) = ,t} , where the sum is over a circular shell defined 

by (n - 1/2) < |k|L/27r < (n + 1/2). 
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FIG. 4: (Color Online) Spatial distribution of the fluid density for the spray-emulsion configuration. The crystal- like ordered 
structure of droplets is evident in fig. (b), where most of droplets are organized into 6- neighbourhood structures. A\ — —4.9 
in all cases, [a]: A 2 = 0.65, nx = ny = 512, t = 500000; [b]: A 2 = 1.15, nx = ny = 512, t = 500000. The spectrum of density 
fluctuations shows a sharp peak corresponding to the typical size of the droplet. For the case (a) we obtain a typical radius of 
R w 8, whereas in the second and more ordered configuration, R « 6.5. 



This first moment of the circularly-averaged structure factor can then be used to assess the characteristic length 
scale of the droplet, R(t) = 2it/k{t), where 

-^-Wti (20) 

Note that k = 1 means R = |r. The right columns in Figs. [2j [4] and [5] show the time evolution of length scales for 
various (Ai, A2). As is well visible from the figures, after 500000 time steps, all configurations have settled down to 
their steady state, except the configuration with A2 ~ A2 C the typical droplet size being a decreasing function of A2. 
It is interesting to notice the growth of macroscopic islands, cutting across the entire computational domain, in the 
emulsion region. This is reflected by a significant build-up of the low-/c region of the spectrum, yet another signature 
of a phase- transit ion behaviour. In figure [6j we show the evolution in time of the typical radius of some configurations. 
The radius is calculated from the circularly averaged structure factor, as described by eq. ([20]) . For the standard 
Shan-Chen case, the radius grows till the maximal value of R ~ L/4, corresponding to a single droplet. In this case, 
the domains grow according to a sub-diffusive power-law R(t) ~ (t — to) a , with a growth exponent a = 1/3 (20j . 
In the other cases, after the transition to the emulsion region, the asymptotic radius attains a much smaller value. 
For A 2 = 0.15, the radius still grows, although more slowly, with a growth exponent a = 1/6, indicating that this 
metastable state will reach the asymptotic single-droplet state in a very long, but finite, time. On the contrary, for 
the other two cases in the emulsion region, the radius does not show any appreciable change over the entire simulation 
time-span. These states appear completely frozen and do not show any visible dynamics towards a more stable state. 
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FIG. 5: (Color online) Spatial distribution of the fluid density for the spray-emulsion configuration. As in previous picture fig. 
[H Ai = —4.9 in all cases and the corresponding Fourier spectra of density fluctuations are reported in the right panel, [a]: 
A 2 = 1.65, nx = ny = 512, t = 500000; [b]: A 2 = 2.15, nx = ny = 512, t = 500000; [c]: A 2 = 2.65, nx = ny = 512, t = 500000. 
Besides the sharp peak centered around the mean size of the droplets, the build-up of a \ow-k component with increasing A 2 is 
well visible, corresponding to the formation of large-scale domains indicating a higher degree of order in the global structure. 
The configuration presented in (c) is strongly reminiscent of a crystal, with very few defects. For these cases, the typical radius 
is estimated as follows: (a) R « 5.1; (b) R « 4.4; (b) R « 3.9. 



V. PSEUDO-POTENTIAL ENERGY EVOLUTION 



The pseudo-potential LB models bears a formal resemblance to dynamic mean-field Ising formulations of magnetic 
systems. Of course, a major difference with respect to Ising systems is that our fluid model is clearly not a Hamiltonian 
one. It is nonetheless of interest to define a pseudo-potential energy, E(t) = E s (t) + E m (t), where 

1 bl 

E s (t) = gZ^'^Z) ( Gl w i + G *P*i) ^si\t) (21) 

x,y 2=0 
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FIG. 6: Time evolution of average domain size R(t) (lattice units) versus time (time steps) for different cases. Curves from top to 
bottom correspond to systems with increased mid-range repulsion ( systems with increased average "surfactant" concentration) . 
The straight line represents a power law R(t) ~ (t — to) 1 / 3 which is typical of diffusive growth. The dashed curve is R(t) ~ 
(t-t ) 1/5 



1 62 

E m(t) = - ^ ^( X > V'i t) ^2 G 2 Pmi l/j(x mi ]t) (22) 



2 

x,y 



are the contributions from the first and second belts, respectively. This definition, suggested by a direct analogy with 
the Ising Hamiltonian H[s] = ^ x Y^ y = x ±i s (y)J( x ^y) s ( x )^ 18 a l so m nne w ^ n the expression of the forces, (pTU|) . 
By expanding ipi in powers of to zeroth order (local-density approximation), we obtain the bulk contribution: 



E bulk = ^-Y,^ 2 ^y) ( 23 ) 



2 



while the next order (weak-gradient approximation) delivers a surface term: 



x,y 

where we used the normalizations in Eqs. (fl2j) and (fl3j) . 

In figure ([7j), the ratio of the global pseudo-energy to the thermal energy E t h = pc 2 s L 2 is shown as a function of time 
for increasing values of the second-belt coupling A2. The figure shows that the steady-state value of the pseudo-energy 
is a monotonically increasing function of A2, the standard SC case (A2 = 0) being the lowest-energy phase-separated 
configuration. The initial rise of the global energy reflects the build-up of surface energy due to interface formation. 
Once such short transient is settled down, the pseudo-energy remains pretty constant in time. Since the "thermal 
energy" E t h is strictly conserved in time the total pseudo-energy, thermal plus potential, may indeed be paralleled 
to a true conserved quantity (Hamiltonian) for most of the time evolution of the system, except a very short initial 
transient. 

In figure ([8]), the time evolution of the ratio of first-belt to second-belt pseudo-energy, for increasing values of 
the parameter A2, is shown. Here again, after a very short transient, the ratio settles down to a constant value, 
which is an increasing function of A2. To be noted that in all cases the ratio is less than 10 percent. Yet, the effect 
on the surface/ volume ratio of the fluid configuration is a very sizeable one, as we shall discuss shortly. The ratio 
between interfacial and bulk components can be estimated as E sur f / E^ik ~ S being the width of interface, 

A the interfacial area and V = L 2 the volume of the simulation box. We have checked that the total volume of 
the liquid phase is the same as in the standard SC case, whereas the interfacial area grows roughly with the scaling 
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FIG. 7: (a): Time evolution of total pseudo-energy E(t) in units of the thermal energy E t h = pc 2 L 2 for increasing value of the 
second belt coupling A2 . The steady-state value of the pseudo-energy is a monotonically increasing function of A2 , the standard 
SC case (A2 = 0) being the lowest-energy phase-separated configuration. In that sense, the standard SC configuration may 
represent a ground state with a discrete spectrum of excited states triggered by increasing discretely A2. (b): Time evolution 
of Etot and Ebuik, Eqs. (20) and (22). The surface energy, given by the difference Etot — E^uik is found to be always positive, as 
it should be. Furthermore, it is possible to see that the surface contribution increases from the Multi-Droplet to the Emulsion 
region (see Fig. 3), consistently with the picture of states which become increasingly excited with increasing A2. In the above 
scales, the two plots in the SC case would almost coincide, since the surface contribution is a mere « 0.001, instead of « 6% 
for the emulsion case. 
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FIG. 8: Time evolution of the ratio Em/E s of the energy associated with the second to first-belt, Eq. (|22|) . The SC line 
(full) represents the ground state. As A2 overcomes A2 C , the energy level increases sensibly due to the interfacial contribution. 
This figure shows the importance of repulsive mid-range interaction in the emulsion region. This interaction is responsible for 
arresting the droplet coalescence sustained by short-range interaction, thereby promoting increased order in the geometrical 
distribution of the droplets. 



relation area/volume ~ n 1 / 2 , n being the number of droplets. This is simply explained in term of mass conservation: 
the volume of a single droplet is given by ttR 2 , while with n droplets, the same volume is given by TrnR^, so that 
R n w Rj y/n. This argument together with the dependence of the number of droplets on the mid-range force G2 shown 
in fig. 2] gives a relation between the final average domain size and the force R ~ G^ 2 . This non-linear dependence 
is seen if fig [6] 

The consistency between theoretical estimation and simulation results has been checked. For instance for the case 
A — 0.15. Typical values are A/V ~ 0.1 and 5/ Ax ~ 10, such that the surface energy should be of the order of 1%. 
This is in line with the actual surface energy, as shown in figure [71 

In figure [9] we show the surface (perimeter in two-dimensions) of the multiphase fluid as a function of time for 
different values of A^. This is seen to go from roughly 5 x 10 -3 of the volume for the SC configuration, up to 0.19 of 
the volume for the emulsion-like configuration, thus showing a factor 40 boost in surface/ volume ratio, even though 
the "potential" energy in the second-shell is just a 10 percent of the energy in the first shell, as shown in fig. [SJ 
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FIG. 9: Total interfacial area, perimeter in 2D, as a function of time for different values of A^. This observable is used 
to monitor the onset of the transition between the multi-droplet region and the emulsion one. In the multi-droplet region 
A2 < A cr , the asymptotic limit is always the equilibrium state with one single droplet, that is the stable minimum in surface 
and, thus, in free-energy. It is clear, however, that the configuration obtained by changing A2 has a different, and yet always 
finite, relaxation time. For A2 > A cr , the area of the liquid phase remains nearly constant in time and the relaxation time 
presents a sharp jump, virtually to an infinite value, thus signalling a phase-transition. The emulsion state is still met ast able, 
but with a life-time much longer than the simulation-time. 

Such a dramatic boost shows that indeed a tiny amount of mid-range repulsion can cause dramatic effects on the 
macroscopic fluid configuration. From this time evolution, it is possible to extract a rough estimate of the equilibrium 
relaxation-time of the system, namely the time necessary to relax to the minimum "free-energy state" (single-droplet). 
For the standard Shan-Chen model, this time has been measured to be t sc ~ 10 3 . 

VI. CONCLUSIONS 

Summarzing, the effects of mid-range repulsion in Lattice Boltzmann models of single-component, non-ideal fluids 
are investigated. The simulations show that mid-range repulsive interactions promote the formation of spray-like, 
multi-droplet configurations, with droplet size directly related to the strength of the repulsive interaction. Our 
results indicate that a small amount of mid-range repulsion can dramatically increase the surface/ volume ratio of the 
multiphase fluid. 

The present approach should offer an useful tool for the computational modelling of complex flow phenomena, 
such as atomization, spray formation and micro-emulsions, break-up phenomena and possibly glassy-like systems as 
well 
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